# Create tibble mapping debrief questions to labels
dd_labels <- tibble(
  var = c("dd2", "dd3", "dd4.1", "dd5", "dd7", "dd8", "dd9"),
  label = c(
    "I said things in order to look good.",
    "I tried to get others to agree with me.",
    "I behaved differently because I was being recorded.",
    "I tried to say what others would agree with.",
    "I tried to influence others.",
    "I was comfortable disagreeing with others.",
    "I let people have other opinions to me."
  )
)

r2_with_dd <- r2_choices_num %>% 
filter(phase == "phase_1") %>% 
mutate(across(matches("^dd\\d(\\.\\d)?$"), ~ if_else(.x == 2, 0, .x))) %>% 
mutate(across(matches("^dd\\d(\\.\\d)?$"), ~ if_else(.x < 0, NA_real_, .x)))
# Wording of dd8, dd9 was reversed after April 8 2023
# dd8 and dd9 only asked after April 8 2023

df %>% count_prop(dd8) # no variation in either one
df %>% count_prop(dd9)

r2_with_dd %>% count_prop(dd7, dd7_label)
# Get dd variables from r2_choices_num
dd_vars <- names(r2_with_dd)[str_detect(names(r2_with_dd), "^dd\\d(\\.\\d)?$")] %>%
  setdiff(c("dd4.1", "dd8", "dd9"))

# Create bar graphs of all dd responses by week
r2_with_dd %>%
  select(survey_date, all_of(dd_vars)) %>%
  pivot_longer(cols = all_of(dd_vars), names_to = "dd_var", values_to = "response") %>%
  left_join(dd_labels, by = c("dd_var" = "var")) %>%
  mutate(
    week = floor_date(survey_date, "week"),
    label = str_glue("{dd_var}: {label}")
  ) %>%
  ggplot(aes(x = week, y = response)) +
  geom_bar(stat = "summary", fun = "mean") +
  facet_wrap(~label, ncol = 2) +
  labs(
    x = "Week",
    y = "Average Agreement", 
    title = "Responses to Discussion Debrief Questions Over Time"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(size = 8)
  )

# Create list of models - one for each dd variable
dd_models <- dd_vars %>% 
  purrr::set_names() %>%
  map(~ feols_custom(
    as.formula(str_glue("r2_choose_trans ~ {.x}")),
    data = r2_with_dd %>% filter(!is.na(r2_choose_trans)),
    fixef = c("stratum_id", "video_type", "delivery_incentive_exp"), 
    cluster = "group_id"
  ))

# Extract heterogeneous effects
het_effects <- dd_models %>%
map_dfr(tidy_90) %>% 
filter(str_detect(term, "^dd\\d(\\.\\d)?$"))


# Calculate proportion of people who agreed (=1) with each debrief question
dd_agreement <- r2_with_dd %>%
  select(all_of(dd_vars)) %>%
  summarise(across(everything(), 
    ~ mean(.x == 1, na.rm = TRUE)
  )) %>%
  pivot_longer(
    everything(),
    names_to = "var",
    values_to = "prop_agree"
  ) %>%
  left_join(dd_labels, by = "var") %>%
  mutate(
    prop_agree_pct = scales::percent(prop_agree, accuracy = 0.1),
    summary = str_glue("{var}: {prop_agree_pct} agreed with '{label}'")
  ) %>%
  arrange(-prop_agree)

# Print results
dd_agreement %>%
  pull(summary) %>%
  cat(sep = "\n")

# Merge dd_labels with het_effects and add agreement proportions to get labels for plot
het_effects_labeled <- het_effects %>%
  inner_join(dd_labels, by = c("term" = "var")) %>%
  # Join with dd_agreement to get proportions
  left_join(dd_agreement %>% select(var, prop_agree_pct), by = c("term" = "var")) %>%
  mutate(
    # Create combined label with question number, text, and agreement proportion
    plot_label = str_glue("{term}: {label} ({prop_agree_pct})")
  )

# Create coefficient plot of heterogeneous effects
ggplot(
  het_effects_labeled,
  aes(x = estimate, y = plot_label)
) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high)
  ) +
  labs(
    x = "Coefficient estimate",
    y = "Debrief question",
    title = "Heterogeneous effects of discussion by debrief responses"
  ) +
  theme_classic()

ggsave("outputs/figs/debrief_het_effects.pdf", width = 8, height = 6)


# First create combined indices for social image and persuasion
r2_with_dd_combined <- r2_with_dd %>%
  mutate(
    # Create social image index (average of dd2 and dd5)
    social_image_index = (dd2 + dd5) / 2,
    # Create persuasion index (average of dd3 and dd7)
    persuasion_index = (dd3 + dd7) / 2
  )

# Create models with combined indices
dd_combined_models <- list(
  "social_image" = feols_custom(
    r2_choose_trans ~ social_image_index,
    data = r2_with_dd_combined %>% filter(!is.na(r2_choose_trans)),
    fixef = c("stratum_id"), 
    cluster = "group_id"
  ),
  "persuasion" = feols_custom(
    r2_choose_trans ~ persuasion_index,
    data = r2_with_dd_combined %>% filter(!is.na(r2_choose_trans)),
    fixef = c("stratum_id"), 
    cluster = "group_id"
  )
)

# Extract effects for plotting
combined_effects <- dd_combined_models %>%
  map_dfr(tidy_90, .id = "mechanism") %>%
  filter(term %in% c("social_image_index", "persuasion_index")) %>%
  mutate(
    mechanism_label = case_when(
      mechanism == "social_image" ~ "Social Image\n(dd2 + dd5)",
      mechanism == "persuasion" ~ "Persuasion\n(dd3 + dd7)"
    )
  )

# Create coefficient plot of combined effects
ggplot(
  combined_effects,
  aes(x = estimate, y = mechanism_label)
) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high)
  ) +
  labs(
    x = "Coefficient estimate",
    y = "Combined mechanism",
    title = "Effects of Combined Discussion Debrief Responses"
  ) +
  theme_classic()

ggsave("outputs/figs/debrief_combined_effects.pdf", width = 8, height = 4)


# Calculate group-level indices (excluding self) for each individual
r2_with_dd_combined <- r2_with_dd_combined %>%
  group_by(group_id) %>%
  mutate(
    # For each mechanism, calculate mean of others in group (excluding self)
    social_image_index_others = (sum(social_image_index) - social_image_index) / (n() - 1),
    persuasion_index_others = (sum(persuasion_index) - persuasion_index) / (n() - 1)
  ) %>%
  ungroup()

# Create models with both individual and group-level indices
dd_combined_models_with_others <- list(
  "social_image" = feols_custom(
    r2_choose_trans ~ social_image_index + social_image_index_others + persuasion_index + persuasion_index_others,
    data = r2_with_dd_combined %>% filter(!is.na(r2_choose_trans)),
    fixef = c("stratum_id"), 
    cluster = "group_id"
  )
)

# Extract effects for plotting
combined_effects_with_others <- dd_combined_models_with_others %>%
  map_dfr(tidy_90, .id = "mechanism") %>%
  filter(term %in% c("social_image_index", "social_image_index_others", 
                     "persuasion_index", "persuasion_index_others")) %>%
  mutate(
    level = if_else(str_detect(term, "_others"), "Group average\n(excluding self)", "Individual"),
    mechanism_clean = if_else(str_detect(term, "social_image"), "Social Image\n(dd2 + dd5)", "Persuasion\n(dd3 + dd7)"),
    plot_label = paste0(mechanism_clean, "\n", level)
  )

# Create coefficient plot showing both individual and group-level effects
ggplot(
  combined_effects_with_others,
  aes(x = estimate, y = plot_label)
) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high)
  ) +
  labs(
    x = "Coefficient estimate",
    y = "Mechanism and level",
    title = "Effects of Individual and Group-Level Discussion Debrief Responses"
  ) +
  theme_classic()

ggsave("outputs/figs/debrief_combined_effects_with_others.pdf", width = 8, height = 6)


# REALLY I SHOULD LOOK AT R1 BEHAVIOUR as the outcome - did they speak up in favour of trans worker or not?

# spoke_in_favour_discuss dataset, includes:
# spoke_pro_trans: 1 if they spoke up in favour of trans worker, 0 if not

df_with_spoke <- df %>% 
left_join(spoke_in_favour_discuss, by = "ind_id", suffix = c("", "_discuss")) %>% 
filter(phase == "phase_1") %>% 
mutate(across(matches("^dd\\d(\\.\\d)?$"), ~ if_else(.x == 2, 0, .x))) %>% 
mutate(across(matches("^dd\\d(\\.\\d)?$"), ~ if_else(.x < 0, NA_real_, .x))) %>% 
filter(!is.na(spoke_pro_trans)) %>% 
mutate(spoke_pro_trans = as.numeric(spoke_pro_trans == 1))

df_with_spoke %>% count_prop(spoke_pro_trans)

# Create models examining difference in dd responses between pro and anti-trans speakers
dd_spoke_diff_models <- dd_vars %>% 
  purrr::set_names() %>%
  map(~ feols_custom(
    as.formula(str_glue("{.x} ~ spoke_pro_trans")),
    data = df_with_spoke,
    fixef = c("stratum_id", "video_type"),
    cluster = "group_id"
  ))

  dd_spoke_diff_models[[1]]$nobs %>% write_stat("outputs/stats/spoke_pro_trans_nobs.tex", digits = 0)

# Extract effects and add labels
dd_spoke_diff_effects <- dd_spoke_diff_models %>%
  map_dfr(tidy_90, .id = "dd_var") %>%
  filter(term == "spoke_pro_trans") %>%
  left_join(dd_labels, by = c("dd_var" = "var")) %>%
  # Join with dd_agreement to get proportions
  left_join(dd_agreement %>% select(var, prop_agree_pct), by = c("dd_var" = "var")) %>%
  mutate(
    plot_label = str_glue("{dd_var}: {label}\n({prop_agree_pct} agreed)")
  )

# Create coefficient plot
ggplot(
  dd_spoke_diff_effects,
  aes(x = estimate, y = plot_label)
) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high)
  ) +
  labs(
    x = "Difference in Response (Spoke Pro-Trans - Did Not)",
    y = "Discussion Debrief Question",
    title = "Differences in Discussion Debrief Responses\nBetween Pro-Trans and Other Speakers"
  ) +
  theme_classic()

ggsave("outputs/figs/debrief_spoke_differences.pdf", width = 8, height = 6)


# Create long format data for plotting
dd_spoke_means <- df_with_spoke %>%
  pivot_longer(
    cols = all_of(dd_vars),
    names_to = "dd_var",
    values_to = "response"
  ) %>%
  # Join with labels
  left_join(dd_labels, by = c("dd_var" = "var")) %>%
  # Calculate means by dd_var and spoke_pro_trans
  mutate(spoke_pro_trans = if_else(spoke_pro_trans >= 0.5, 1, 0)) %>%
  group_by(dd_var, label, spoke_pro_trans) %>%
  summarise(
    mean_response = mean(response, na.rm = TRUE),
    se = sd(response, na.rm = TRUE) / sqrt(n()),
    .groups = "drop"
  ) %>%
  mutate(
    spoke_label = if_else(spoke_pro_trans >= 0.5, 
                         "Spoke up for trans workers", 
                         "Did not speak up for trans workers")
  )



# --- Prepare data for improved bar plot ---

# Only include debrief questions in the two mechanisms:
plot_vars <- c("dd2", "dd3", "dd5", "dd7")

# Start with the raw computed means for each question and speaking group:
dd_spoke_plot <- dd_spoke_means %>%
  filter(dd_var %in% plot_vars) %>%
  mutate(
    mechanism = case_when(
      dd_var %in% c("dd2", "dd5") ~ "Social image",
      dd_var %in% c("dd3", "dd7") ~ "Persuasion"
    ),
    # Create a facet label that combines the dd var and the full question text
    facet_label = str_glue("{label}")
  ) %>%
  select(dd_var, mechanism, facet_label, spoke_label, mean_response)

# Extract baseline values (raw means for those who did NOT speak up):
baseline_df <- dd_spoke_plot %>%
  filter(spoke_label == "Did not speak up for trans workers") %>%
  rename(baseline_mean = mean_response) %>%
  select(dd_var, mechanism, facet_label, baseline_mean)

# Extract regression estimates for each dd question (remember this is a model
# run with feols_custom() and coefficients extracted using tidy_90(), which we assume
# now returns 95% CIs)
dd_reg <- dd_spoke_diff_effects %>% 
  filter(dd_var %in% plot_vars) %>%
  select(dd_var, estimate, conf.low, conf.high, p.value) %>%
  mutate(
    stars = case_when(
      p.value < 0.01 ~ "***",
      p.value < 0.05 ~ "**",
      p.value < 0.1  ~ "*",
      TRUE ~ ""
    )
  )

# Combine the regression estimates with the baseline so that the "Spoke up for trans workers"
# mean is computed as: raw mean (for non-speakers) + coefficient.
reg_adjusted <- baseline_df %>%
  left_join(dd_reg, by = "dd_var") %>%
  mutate(
    mean_spoke = baseline_mean + estimate,
    ci_low = baseline_mean + conf.low,
    ci_high = baseline_mean + conf.high
  ) %>%
  select(dd_var, mechanism, facet_label, mean_spoke, ci_low, ci_high, stars, p.value)

# Now create the final plotting dataset with one row for each bar:
baseline_bars <- baseline_df %>%
  mutate(speaker = "Did not speak up for trans workers",
         mean_value = baseline_mean) %>%
  select(dd_var, mechanism, facet_label, speaker, mean_value)

spoke_bars <- reg_adjusted %>%
  mutate(speaker = "Spoke up for trans workers",
         mean_value = mean_spoke) %>%
  select(dd_var, mechanism, facet_label, speaker, mean_value, ci_low, ci_high, stars, p.value)

plot_df <- bind_rows(baseline_bars, spoke_bars) %>%
  mutate(speaker = factor(speaker, levels = c("Did not speak up for trans workers",
                                                "Spoke up for trans workers")))

# Compute positions for significance brackets.
# For each dd_var (i.e. per facet) extract the top y-value from the "Spoke up" bar and add an offset.
signif_df <- plot_df %>%
  filter(speaker == "Spoke up for trans workers") %>%
  group_by(dd_var, mechanism, facet_label, speaker) %>%
  summarize(
    top_y = max(ci_high, mean_value, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  left_join(dd_reg, by = "dd_var") %>%
  mutate(
    # Position the bracket a little above the top of the error bar/mean:
    y.position = top_y + 0.05,
    # Create a label with the p-value and significance stars
    label = str_glue("p={formatC(p.value, format = 'f', digits = 2)} {stars}") %>% as.character()
  )

# First update the facet labels to combine mechanism and wrapped text
plot_df <- plot_df %>%
  mutate(
    # Create combined label with mechanism and wrapped text
    facet_label = str_wrap(str_glue("{mechanism}:\n{facet_label}"), width = 35)
  )

# Update signif_df to match the new labels
signif_df <- signif_df %>%
  mutate(
    facet_label = str_wrap(str_glue("{mechanism}:\n{facet_label}"), width = 35)
  )

# Update speaker labels
plot_df <- plot_df %>%
  mutate(
    speaker = factor(
      speaker,
      levels = c("Did not speak up for trans workers", "Spoke up for trans workers"),
      labels = c("Did not speak up\nin favor of trans", "Spoke up\nin favor of trans")
    )
  )

  

# Force the ordering of facets so that Social image panels appear on top and Influence panels below
plot_df <- plot_df %>%
  mutate(
    mechanism = factor(mechanism, levels = c("Persuasion", "Social image"))
  ) %>%
  arrange(mechanism) %>%
  mutate(
    facet_label = factor(as.character(facet_label), levels = c(
      "Social image: I tried to say what\nothers would agree with.",
      "Social image: I said things in\norder to look good.",
      "Persuasion: I tried to influence\nothers.",
      "Persuasion: I tried to get others\nto agree with me." 
    )) %>% 
    fct_rev()
  ) %>% 
  arrange(facet_label)


ggplot(plot_df, aes(x = speaker, y = mean_value, fill = speaker)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), alpha = 0.8, colour = "darkgrey") +
  # Add value labels at fixed y position
  geom_text(aes(label = sprintf("%.2f", mean_value), y = 0.1), 
            size = 3, color = "black") +
  # Only add error bars for "Spoke up for trans workers" (the regression estimates):
  geom_errorbar(data = filter(plot_df, speaker == "Spoke up\nin favor of trans"),
                aes(ymin = ci_low, ymax = ci_high),
                width = 0.2, position = position_dodge(width = 0.9)) +
  scale_fill_manual(values = c("Did not speak up\nin favor of trans" = "#E0E0E0", 
                              "Spoke up\nin favor of trans" = "#F2A365")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  facet_wrap(~ facet_label, ncol = 2) +
  labs(
    x = NULL,
    y = "Mean Agreement"
  ) +
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 0, hjust = 0.5),
    panel.grid = element_blank()
  ) +
  # --- Add slightly narrower significance brackets above each pair of bars ---
  geom_segment(data = signif_df, 
               aes(x = 0.8, xend = 2.2, y = y.position, yend = y.position)) +
  geom_segment(data = signif_df, 
               aes(x = 0.8, xend = 0.8, y = y.position - 0.03, yend = y.position)) +
  geom_segment(data = signif_df, 
               aes(x = 2.2, xend = 2.2, y = y.position - 0.03, yend = y.position)) +
  geom_text(data = signif_df, 
            aes(x = 1.5, y = y.position + 0.02, label = label),
            size = 3, vjust = 0)

ggsave("outputs/figs/debrief_spoke_differences_bar.pdf", width = 5, height = 6)


# Export mean of spoke_pro_trans as a stat
df_with_spoke %>% get_mean("spoke_pro_trans") %>% write_percentage("outputs/stats/spoke_pro_trans_mean.tex")


dd_spoke_diff_models[["dd3"]] %>% get_p_val("spoke_pro_trans") %>% write_p_val("outputs/stats/spoke_pro_trans_dd3_pval.tex")
dd_spoke_diff_models[["dd3"]] %>% get_coeff("spoke_pro_trans") %>% times_100 %>% write_stat("outputs/stats/spoke_pro_trans_dd3_coeff.tex", digits = 1)
((dd_spoke_diff_models[["dd3"]] %>% get_coeff("spoke_pro_trans")) / (dd_spoke_means %>% filter(dd_var == "dd3", spoke_pro_trans == 0) %>% pull(mean_response))) %>% write_percentage("outputs/stats/spoke_pro_trans_dd3_prop_increase.tex")


# Export p-values of other dd variables
dd_spoke_diff_models[["dd2"]] %>% get_p_val("spoke_pro_trans") %>% write_p_val("outputs/stats/spoke_pro_trans_dd2_pval.tex")
dd_spoke_diff_models[["dd5"]] %>% get_p_val("spoke_pro_trans") %>% write_p_val("outputs/stats/spoke_pro_trans_dd5_pval.tex")
dd_spoke_diff_models[["dd7"]] %>% get_p_val("spoke_pro_trans") %>% write_p_val("outputs/stats/spoke_pro_trans_dd7_pval.tex")
